library(lattice)
library(foreign)
library(Zelig)

cat('\nEstimating Senate election effects ...\n')

#
#Read in the state partisan composition data
#
wright<- read.dta("exitpolls/cbs7603_percents.dta",convert.underscore=F)

#
# Make data contest data table for use in producing counterfactual first differences
#
senate.contests <- sort(unique(aeh.senate$contest))
n.senate.contests <- length(senate.contests)
senate.contests.prediction.data <- matrix(NA,nrow=n.senate.contests,ncol=10)
colnames(senate.contests.prediction.data) <- c("id.contest","face.chal","face.inc","cook","inc.yrsinoffice","inc.age","chalparty.pid.pct","incparty.pid.pct","inc.margin","contest.label")
for (i in 1:n.senate.contests){
  if (all(c(0,1) %in% aeh.senate$incumbent[aeh.senate$contest==senate.contests[i]]) & senate.contests[i] %in% aeh.senate$contest){
    senate.contests.prediction.data[i,1] <- senate.contests[i]
    senate.contests.prediction.data[i,2] <- aeh.senate$comp[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==0]
    senate.contests.prediction.data[i,3] <- aeh.senate$comp[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]
    cook.label <- aeh.senate$cook[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]
    senate.contests.prediction.data[i,4] <- ifelse(toupper(cook.label)=="TOSSUPDEM" | toupper(cook.label)=="TOSSUPREP",3,ifelse(toupper(cook.label)=="LEANDEM" | toupper(cook.label)=="LEANREP",2,ifelse(toupper(cook.label)=="LIKELYDEM" | toupper(cook.label)=="LIKELYREP",1,ifelse(toupper(cook.label)=="SOLIDDEM" | toupper(cook.label)=="SOLIDREP",0,NA))))
    senate.contests.prediction.data[i,5] <- aeh.senate$year[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]-aeh.senate$IncFirstYrInOffice[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]
    senate.contests.prediction.data[i,6] <- aeh.senate$year[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]-aeh.senate$IncBirthYr[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]
    inc.is.dem <- ifelse(aeh.senate$party[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]=="dem",1,ifelse(aeh.senate$party[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]=="rep",0,NA))
    wright.dem <- ifelse(aeh.senate$year[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]<=2003,
                    wright$wtd_pty3[wright$year==aeh.senate$year[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1] & substr(aeh.senate$office[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1],5,6)==wright$stateid],
                    wright$wtd_pty3[wright$year==2003 & substr(aeh.senate$year[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1],5,6)==wright$stateid])
    wright.rep <- ifelse(aeh.senate$year[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]<=2003,
                    wright$wtd_pty1[wright$year==aeh.senate$year[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1] & substr(aeh.senate$office[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1],5,6)==wright$stateid],
                    wright$wtd_pty1[wright$year==2003 & substr(aeh.senate$year[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1],5,6)==wright$stateid])
    senate.contests.prediction.data[i,7] <- ifelse(inc.is.dem==1,wright.rep,wright$wtd_pty3)
    senate.contests.prediction.data[i,8] <- ifelse(inc.is.dem==1,wright.dem,wright.rep)
    total.votes <- aeh.senate$votes[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]+aeh.senate$votes[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==0]
    senate.contests.prediction.data[i,9] <- (aeh.senate$votes[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]/total.votes) - (aeh.senate$votes[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==0]/total.votes)
    senate.contests.prediction.data[i,10] <- aeh.senate$office[aeh.senate$contest==senate.contests[i] & aeh.senate$incumbent==1]
  }
}
senate.contests.prediction.data <- as.data.frame(na.omit(senate.contests.prediction.data),stringsAsFactors=F)

senate.contests.prediction.data[1,]

#
# Zelig predicted values analysis
#

# Since we only care about the point estimate, no reason to cluster.
z.out <- zelig(as.formula(senate.formula), model="probit",
                data=exit.senate)
face.inc <- median(c(senate.dem$comp,senate.rep$comp),na.rm=T)
inc.yrsinoffice <- mean(exit.senate$inc.yrsinoffice)
inc.age <- mean(exit.senate$inc.age)

first.dif.sim <- function(chal.copart=0,inc.copart=0,face.chal,face.inc,cook,
                        inc.yrsinoffice,inc.age){
  x.out <- setx(z.out)
  x.out$cook <- cook
  x.out$inc.yrsinoffice <- inc.yrsinoffice
  x.out$inc.age <- inc.age
  x.out$shares.chal.party <- chal.copart
  x.out$shares.inc.party <- inc.copart
  x.out$face.inc  <- face.inc
  x.out_median <- x.out
  x.out_median$face.chal <-  quantile(c(senate.dem$comp,senate.rep$comp),.50,na.rm=T)
  x.out_actual <- x.out
  x.out_actual$face.chal <- face.chal
  s.out <- sim(z.out, x = x.out_median, x1 = x.out_actual)
  return(s.out$qi$fd)
}

n.contests <-  nrow(senate.contests.prediction.data)
effect.out <- matrix(NA,nrow=n.contests,ncol=3)
for (k in 1:n.contests){
  cat("Contest",k,"\n") ; flush.console()
  independent <- first.dif.sim(chal.copart=0,inc.copart=0,face.chal=as.numeric(senate.contests.prediction.data$face.chal[k]),face.inc=as.numeric(senate.contests.prediction.data$face.inc[k]),cook=as.numeric(senate.contests.prediction.data$cook[k]),inc.yrsinoffice=as.numeric(senate.contests.prediction.data$inc.yrsinoffice[k]),inc.age=as.numeric(senate.contests.prediction.data$inc.age[k]))
  chal <- first.dif.sim(chal.copart=1,inc.copart=0,face.chal=as.numeric(senate.contests.prediction.data$face.chal[k]),face.inc=as.numeric(senate.contests.prediction.data$face.inc[k]),cook=as.numeric(senate.contests.prediction.data$cook[k]),inc.yrsinoffice=as.numeric(senate.contests.prediction.data$inc.yrsinoffice[k]),inc.age=as.numeric(senate.contests.prediction.data$inc.age[k]))
  inc <- first.dif.sim(chal.copart=0,inc.copart=1,face.chal=as.numeric(senate.contests.prediction.data$face.chal[k]),face.inc=as.numeric(senate.contests.prediction.data$face.inc[k]),cook=as.numeric(senate.contests.prediction.data$cook[k]),inc.yrsinoffice=as.numeric(senate.contests.prediction.data$inc.yrsinoffice[k]),inc.age=as.numeric(senate.contests.prediction.data$inc.age[k]))
  chal.pct <- as.numeric(senate.contests.prediction.data$chalparty.pid.pct[k])/100
  inc.pct <- as.numeric(senate.contests.prediction.data$incparty.pid.pct[k])/100
  ind.pct <- 1-chal.pct-inc.pct
  effect <- chal.pct*chal+inc.pct*inc+ind.pct*independent
  effect.out[k,] <- c(mean(effect),quantile(effect,c(.025,.975)))
}

effects.data <- list()
effects.data$chal.effect <- -1*effect.out[,1] #score margin in chal direction
effects.data$contest <- paste(substr(senate.contests.prediction.data$contest.label,5,6),
                            substr(senate.contests.prediction.data$contest.label,3,4))
effects.data$inc.margin <- as.numeric(senate.contests.prediction.data$inc.margin)
effects.data <- as.data.frame(effects.data,stringsAsFactors=F)

# Identify effect sizes greater than margin
changes.outcome <- abs(effects.data$chal.effect)>abs(effects.data$inc.margin) & !is.na(effects.data$inc.margin)
cat("Elections that would have had a different outcome:\n")
print(effects.data[changes.outcome,])
cat("Santorum's face effect:\n")
print(effects.data[effects.data$contest=='PA 94',])

#
#Data for the NM2000 example used in the text.
#
independent <- first.dif.sim(chal.copart=0,inc.copart=0,face.chal=as.numeric(senate.contests.prediction.data$face.chal[64]),face.inc=as.numeric(senate.contests.prediction.data$face.inc[64]),cook=as.numeric(senate.contests.prediction.data$cook[64]),inc.yrsinoffice=as.numeric(senate.contests.prediction.data$inc.yrsinoffice[64]),inc.age=as.numeric(senate.contests.prediction.data$inc.age[64]))
chal <- first.dif.sim(chal.copart=1,inc.copart=0,face.chal=as.numeric(senate.contests.prediction.data$face.chal[64]),face.inc=as.numeric(senate.contests.prediction.data$face.inc[64]),cook=as.numeric(senate.contests.prediction.data$cook[64]),inc.yrsinoffice=as.numeric(senate.contests.prediction.data$inc.yrsinoffice[64]),inc.age=as.numeric(senate.contests.prediction.data$inc.age[64]))
inc <- first.dif.sim(chal.copart=0,inc.copart=1,face.chal=as.numeric(senate.contests.prediction.data$face.chal[64]),face.inc=as.numeric(senate.contests.prediction.data$face.inc[64]),cook=as.numeric(senate.contests.prediction.data$cook[64]),inc.yrsinoffice=as.numeric(senate.contests.prediction.data$inc.yrsinoffice[64]),inc.age=as.numeric(senate.contests.prediction.data$inc.age[64]))
chal.pct <- as.numeric(senate.contests.prediction.data$chalparty.pid.pct[64])/100
inc.pct <- as.numeric(senate.contests.prediction.data$incparty.pid.pct[64])/100
ind.pct <- 1-chal.pct-inc.pct
effect <- chal.pct*chal+inc.pct*inc+ind.pct*independent
cat("Challenger face chal effect in NM2000 is: ",-1*round(mean(effect),4),"\n",
    "For independents, challenger face chal effect in NM2000 is: ",-1*round(mean(independent),4),"\n",
    "For incumbent copartisans (Democrats), challenger face chal effect in NM2000 is: ",-1*round(mean(inc),4),"\n",
    "For challenger copartisans (Republicans), challenger face chal effect in NM2000 is: ",-1*round(mean(chal),4),"\n",
    "In 2000, the percentage of people in NM who aligned with the challenger's party was: ",round(chal.pct,4),"\n",
    "In 2000, the percentage of people in NM who aligned with the incumbent's party was: ",round(inc.pct,4),"\n",
    "In 2000, the percentage of people in NM who aligned as independents was: ",round(ind.pct,4),"\n")
    
   

if(toEPS) {postscript("TableAndFiguresOutput/ChallengerEffects.eps",
                        width=6,height=6,horizontal=F)}
par(mar=c(4.1,1.1,.1,1.1))
max.effect <- round(max(abs(effects.data$chal.effect)),2)
plot(effects.data$chal.effect,1:nrow(effects.data),
    type="n",axes=F,xlim=c(-max.effect,max.effect),ann=F)
effects.data = effects.data[order(effects.data$chal.effect),]
offset = .005
for(i in 1:nrow(effects.data)){
    if(effects.data$chal.effect[i]>0){
        lines(x=c(offset,effects.data$chal.effect[i]+offset),y=c(i,i),lwd=2)
    } else {
        lines(x=c(-offset,effects.data$chal.effect[i]-offset),y=c(i,i),lwd=2)
    }
}
posticks = seq(from=offset,to=max.effect,by=offset*2)
negticks = seq(from=-max.effect+offset,to=-offset,by=offset*2)
neglabs = as.character(100*(negticks+offset))
neglabs[length(neglabs)]="0"
axis(side=1,at=c(negticks,posticks),cex.axis=.8,
    labels=as.character(c(neglabs,100*(posticks-offset))) )
text(x=rep(0,nrow(effects.data)),y=seq(1:nrow(effects.data)),
    labels = effects.data$contest,cex=.4)
mtext(text="Effect on Challenger Vote",at = 0,side=1,line = 3)
if (toEPS) dev.off()
    
